Random integer (discrete uniform) distribution (randint)

The randint distribution models an integer drawn uniformly at random from a finite range of consecutive integers.

In SciPy, it corresponds to scipy.stats.randint(low, high) with support: $\( \{\texttt{low},\ \texttt{low}+1,\ \ldots,\ \texttt{high}-1\}. \)$

Learning goals

• Recognize when a bounded integer uniform model is appropriate.

• Write the PMF and CDF carefully (including the half-open interval convention).

• Compute mean/variance/skewness/kurtosis and entropy.

• Derive the likelihood and the MLE for the bounds.

• Implement sampling from scratch (NumPy-only) and validate it by simulation.

• Use scipy.stats.randint for PMF/CDF/sampling and understand how to fit with scipy.stats.fit.

Prerequisites

• Basic probability (PMF/CDF), expectation, and variance

• Familiarity with sums like \(\sum_{k=0}^{n-1} k\) and \(\sum_{k=0}^{n-1} k^2\)

• Comfort with logs and basic optimization intuition

Notebook roadmap

1. Title & Classification

2. Intuition & Motivation

3. Formal Definition

4. Moments & Properties

5. Parameter Interpretation

6. Derivations (Expectation, Variance, Likelihood)

7. Sampling & Simulation (NumPy-only)

8. Visualization (PMF, CDF, Monte Carlo)

9. SciPy Integration

10. Statistical Use Cases

11. Pitfalls

12. Summary

import numpy as np

import plotly.graph_objects as go
import os
import plotly.io as pio

import scipy
from scipy import stats
from scipy.special import logsumexp

pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")

np.set_printoptions(precision=4, suppress=True)

SEED = 7
rng = np.random.default_rng(SEED)

import sys
import plotly

print("Python:", sys.version.split()[0])
print("NumPy:", np.__version__)
print("SciPy:", scipy.__version__)
print("Plotly:", plotly.__version__)
print("Seed:", SEED)

Python: 3.12.9
NumPy: 1.26.2
SciPy: 1.15.0
Plotly: 6.5.2
Seed: 7

1) Title & Classification

• Name: randint (random integer / discrete uniform on a contiguous integer interval)

• Type: Discrete

• Support (SciPy convention): $\(x \in \{\ell,\ \ell+1,\ \ldots,\ h-1\}\)\( where \)\ell = \texttt{low}\( and \)h = \texttt{high}$.

• Parameter space: $\(\ell \in \mathbb{Z},\quad h \in \mathbb{Z},\quad h > \ell.\)$

It is often convenient to define the number of possible outcomes: $\( n = h - \ell. \)\( Then the support contains exactly \)n$ integers.

Note on SciPy’s loc: many SciPy distributions also accept a loc shift. For randint, loc shifts the support by addition: $\( X \sim \texttt{randint}(\ell, h, \texttt{loc})\quad\Rightarrow\quad \text{support } \{\ell+\texttt{loc}, \ldots, h-1+\texttt{loc}\}. \)$ In this notebook we focus on the common loc=0 case unless stated otherwise.

2) Intuition & Motivation

This is the discrete analogue of the continuous uniform distribution: all allowed integer values are equally likely.

• What it models: an integer outcome with no reason to prefer one value over another within known bounds.

  • Examples: a fair die roll (with an appropriate mapping), choosing a random index into an array, random day-of-week encoding, random A/B group assignment when the groups are equally sized by design.

• Typical real-world use cases:

  • Random indexing / shuffling: picking a random element from a list.

  • Simulation: drawing random IDs, discrete time steps, or random augmentation choices.

  • Bounded non-informative priors in Bayesian modeling for integer-valued parameters (e.g., a model index among \(\{1,\dots,M\}\)).

Relations to other distributions

• Continuous uniform: if \(U \sim \mathrm{Uniform}(0,1)\) and you set \(X = \ell + \lfloor nU \rfloor\), then \(X\) is randint on \(\{\ell,\dots,h-1\}\).

• Categorical: randint is a categorical distribution over consecutive integers with equal probabilities.

• Bernoulli: when \(n=2\) (two consecutive integers), this is a Bernoulli distribution after a simple re-labeling.

• Discrete uniform on \(\{1,\dots,N\}\): setting low=1, high=N+1 gives the familiar “uniform from 1 to N” model.

3) Formal Definition

Let \(X\) be uniformly distributed on the integers \(\{\ell, \ell+1, \ldots, h-1\}\) where \(\ell,h\in\mathbb{Z}\) and \(h>\ell\). Let \(n = h-\ell\).

PMF

For integer \(k\), $\( \mathbb{P}(X=k) = \begin{cases} \frac{1}{n} & k \in \{\ell,\ell+1,\ldots,h-1\}\\ 0 & \text{otherwise.} \end{cases} \)$

CDF

Because this is a discrete distribution, the CDF is a step function. For real \(x\), $\( F(x)=\mathbb{P}(X\le x)= \begin{cases} 0 & x < \ell \\ \frac{\lfloor x\rfloor-\ell+1}{n} & \ell \le x < h-1 \\ 1 & x \ge h-1. \end{cases} \)$

An equivalent “clipped” form (useful for implementation) is: $\( F(x)=\mathrm{clip}\left(\frac{\lfloor x\rfloor-\ell+1}{n},\ 0,\ 1\right). \)$

def _as_int_like(name: str, value) -> int:
    if isinstance(value, (int, np.integer)):
        return int(value)
    if isinstance(value, (float, np.floating)) and float(value).is_integer():
        return int(value)
    raise TypeError(f"{name} must be an integer (or integer-valued float), got {value!r}")


def validate_low_high(low, high) -> tuple[int, int]:
    low = _as_int_like("low", low)
    high = _as_int_like("high", high)
    if high <= low:
        raise ValueError(f"Require high > low, got low={low}, high={high}")
    return low, high


def randint_support(low, high) -> np.ndarray:
    low, high = validate_low_high(low, high)
    return np.arange(low, high)


def randint_pmf(x, low, high):
    """PMF of randint(low, high) evaluated at x."""
    low, high = validate_low_high(low, high)
    n = high - low

    x = np.asarray(x)
    is_int = x == np.floor(x)
    in_support = is_int & (x >= low) & (x < high)

    return np.where(in_support, 1.0 / n, 0.0).astype(float)


def randint_cdf(x, low, high):
    """CDF of randint(low, high) evaluated at x (step function)."""
    low, high = validate_low_high(low, high)
    n = high - low

    x = np.asarray(x, dtype=float)
    m = np.floor(x)
    cdf = (m - low + 1.0) / n

    return np.clip(cdf, 0.0, 1.0)


low_demo, high_demo = 2, 7
x_demo = np.array([1, 2, 3, 6, 7, 8], dtype=float)
print("support:", randint_support(low_demo, high_demo))
print("x:", x_demo)
print("pmf:", randint_pmf(x_demo, low_demo, high_demo))
print("cdf:", randint_cdf(x_demo, low_demo, high_demo))

support: [2 3 4 5 6]
x: [1. 2. 3. 6. 7. 8.]
pmf: [0.  0.2 0.2 0.2 0.  0. ]
cdf: [0.  0.2 0.4 1.  1.  1. ]

4) Moments & Properties

Let \(n=h-\ell\).

Mean and variance

Because the distribution is uniform on a finite set, all moments exist.

• Mean: $\(\mathbb{E}[X] = \frac{\ell + (h-1)}{2} = \frac{\ell + h - 1}{2}.\)$

• Variance: $\(\mathrm{Var}(X) = \frac{n^2-1}{12}.\)$

Skewness and kurtosis

• Skewness: \(0\) (the distribution is symmetric about its mean).

• Excess kurtosis (Fisher kurtosis): $\(\gamma_2 = -\frac{6(n^2+1)}{5(n^2-1)}.\)$

MGF and characteristic function

For \(t\ne 0\), the moment generating function is a finite geometric series: $\( M_X(t)=\mathbb{E}[e^{tX}] = \frac{1}{n}\sum_{k=\ell}^{h-1} e^{tk} = \frac{e^{t\ell}\,(1-e^{tn})}{n(1-e^{t})}. \)\( (And \)M_X(0)=1$ by continuity.)

The characteristic function is the same expression with \(t=i\omega\): $\( \varphi_X(\omega) = \mathbb{E}[e^{i\omega X}] = \frac{e^{i\omega\ell}(1-e^{i\omega n})}{n(1-e^{i\omega})} = e^{i\omega(\ell + (n-1)/2)}\,\frac{\sin(n\omega/2)}{n\sin(\omega/2)}. \)$

Entropy

Because the distribution is uniform over \(n\) outcomes, $\( H(X) = -\sum_k \mathbb{P}(X=k)\log \mathbb{P}(X=k) = \log n \)\( (in **nats**; use \)\log_2$ for bits).

def randint_moments(low, high):
    low, high = validate_low_high(low, high)
    n = high - low

    mean = 0.5 * (low + high - 1)
    var = (n**2 - 1) / 12
    skew = 0.0 if n > 1 else np.nan
    excess_kurtosis = (-6 * (n**2 + 1) / (5 * (n**2 - 1))) if n > 1 else np.nan
    entropy = float(np.log(n))

    return mean, var, skew, excess_kurtosis, entropy


def randint_mgf(t, low, high):
    """MGF M_X(t) for randint(low, high) using a stable expm1 ratio."""
    low, high = validate_low_high(low, high)
    n = high - low

    t = np.asarray(t, dtype=float)
    out = np.empty_like(t, dtype=float)

    mask0 = t == 0
    out[mask0] = 1.0

    tt = t[~mask0]
    out[~mask0] = np.exp(tt * low) * (np.expm1(tt * n) / np.expm1(tt)) / n

    return out


low, high = 2, 7
mean, var, skew, ex_kurt, ent = randint_moments(low, high)
rv = stats.randint(low, high)

print("Formulas:")
print("  mean:", mean)
print("  var :", var)
print("  skew:", skew)
print("  excess kurtosis:", ex_kurt)
print("  entropy (nats):", ent)

print()
print("SciPy:")
print("  stats(mvsk):", rv.stats(moments="mvsk"))
print("  entropy     :", rv.entropy())

# Numerical check of the MGF against Monte Carlo
n_mc = 200_000
t = 0.3
x_mc = rv.rvs(size=n_mc, random_state=rng)
mgf_emp = np.mean(np.exp(t * x_mc))
mgf_theo = float(randint_mgf(t, low, high))

print()
print("MGF check at t=0.3:")
print("  empirical E[exp(tX)]:", mgf_emp)
print("  theoretical M_X(t)  :", mgf_theo)

Formulas:
  mean: 4.0
  var : 2.0
  skew: 0.0
  excess kurtosis: -1.3
  entropy (nats): 1.6094379124341003

SciPy:
  stats(mvsk): (4.0, 2.0, 0.0, -1.3)
  entropy     : 1.6094379124341003

MGF check at t=0.3:
  empirical E[exp(tX)]: 3.6290750483186573
  theoretical M_X(t)  : 3.626635073807003

5) Parameter Interpretation

randint has two boundary parameters:

• \(\ell=\texttt{low}\) sets the left endpoint (included).

• \(h=\texttt{high}\) sets the right endpoint (excluded).

• \(n=h-\ell\) is the support size (how many integers are possible).

How parameters change the distribution:

• Changing \(\ell\) with fixed \(n\) shifts the distribution left/right but does not change its shape.

• Increasing \(n\) (moving \(h\) farther from \(\ell\)) spreads mass over more integers, decreasing each point probability from \(1/n\).

A useful mental model is: $\( X = \ell + Y,\qquad Y \sim \text{Uniform on } \{0,1,\ldots,n-1\}. \)\( Shifts affect the mean; the variance depends only on \)n$.

6) Derivations

Let \(X\) be uniform on \(\{\ell,\ell+1,\ldots,h-1\}\) and let \(n=h-\ell\).

Expectation

Write \(X=\ell+Y\) where \(Y\) is uniform on \(\{0,1,\ldots,n-1\}\). Then $\( \mathbb{E}[X]=\ell+\mathbb{E}[Y]. \)\( Compute \)\( \mathbb{E}[Y]=\frac{1}{n}\sum_{k=0}^{n-1} k = \frac{1}{n}\cdot\frac{(n-1)n}{2}=\frac{n-1}{2}, \)\( so \)\( \mathbb{E}[X]=\ell+\frac{n-1}{2}=\frac{\ell+h-1}{2}. \)$

Variance

Because variance is shift-invariant, \(\mathrm{Var}(X)=\mathrm{Var}(Y)\). Use \(\mathrm{Var}(Y)=\mathbb{E}[Y^2]-(\mathbb{E}[Y])^2\).

First, $\( \mathbb{E}[Y^2]=\frac{1}{n}\sum_{k=0}^{n-1} k^2 =\frac{1}{n}\cdot\frac{(n-1)n(2n-1)}{6} = \frac{(n-1)(2n-1)}{6}. \)\( Then \)\( \mathrm{Var}(Y)=\frac{(n-1)(2n-1)}{6}-\left(\frac{n-1}{2}\right)^2 =\frac{n^2-1}{12}. \)$

Likelihood (i.i.d. sample)

Let \(x_1,\ldots,x_m\) be i.i.d. draws from randint(\ell,h). The likelihood is $\( L(\ell,h; x_{1:m})=\prod_{i=1}^m \mathbb{P}(X=x_i) =\begin{cases} \left(\frac{1}{h-\ell}\right)^m & \text{if all } x_i \in \{\ell,\ldots,h-1\} \\ 0 & \text{otherwise.} \end{cases} \)$

So the log-likelihood (when feasible) is $\( \ell(\ell,h) = -m\log(h-\ell). \)\( To maximize it, we want the **smallest interval** that still contains the data. Let \)x_{\min}=\min_i x_i\( and \)x_{\max}=\max_i x_i\(. The unique MLE under the half-open support convention is: \)\( \hat{\ell}=x_{\min},\qquad \hat{h}=x_{\max}+1. \)$

def randint_log_likelihood(low, high, x) -> float:
    """Log-likelihood for i.i.d. data x under randint(low, high)."""
    low, high = validate_low_high(low, high)
    x = np.asarray(x)
    if x.size == 0:
        raise ValueError("x must be non-empty")

    if not np.all(x == np.floor(x)):
        raise ValueError("x must contain integer-valued observations")

    x_min = x.min()
    x_max = x.max()
    if (x_min < low) or (x_max >= high):
        return -np.inf

    n = high - low
    return -x.size * float(np.log(n))


def randint_mle(x) -> tuple[int, int]:
    x = np.asarray(x)
    if x.size == 0:
        raise ValueError("x must be non-empty")
    if not np.all(x == np.floor(x)):
        raise ValueError("x must contain integer-valued observations")

    low_hat = int(x.min())
    high_hat = int(x.max()) + 1
    return low_hat, high_hat


# Simulate data and visualize how the likelihood prefers the tightest interval.
low_true, high_true = 3, 11
m = 60
x = stats.randint.rvs(low_true, high_true, size=m, random_state=rng)

low_hat, high_hat = randint_mle(x)
print("true (low, high):", (low_true, high_true))
print("MLE  (low, high):", (low_hat, high_hat))

low_grid = np.arange(low_hat - 6, low_hat + 1)
high_grid = np.arange(high_hat, high_hat + 7)

ll = np.empty((low_grid.size, high_grid.size), dtype=float)
for i, lo in enumerate(low_grid):
    for j, hi in enumerate(high_grid):
        ll[i, j] = randint_log_likelihood(lo, hi, x)

ll_plot = ll.copy()
ll_plot[~np.isfinite(ll_plot)] = np.nan

fig = go.Figure(
    data=go.Heatmap(
        z=ll_plot,
        x=high_grid,
        y=low_grid,
        colorbar_title="log L",
    )
)
fig.add_trace(
    go.Scatter(
        x=[high_hat],
        y=[low_hat],
        mode="markers",
        marker=dict(color="red", size=10),
        name="MLE",
    )
)
fig.update_layout(
    title="Log-likelihood over candidate (low, high) intervals",
    xaxis_title="high (exclusive)",
    yaxis_title="low",
)
fig.show()

true (low, high): (3, 11)
MLE  (low, high): (3, 11)

7) Sampling & Simulation (NumPy-only)

A simple sampler uses a uniform random variable and a floor operation.

1. Draw \(U \sim \mathrm{Uniform}(0,1)\).

2. Return $\(X = \ell + \lfloor nU \rfloor,\quad n=h-\ell.\)$

Why this works: \(\lfloor nU \rfloor\) takes values in \(\{0,\ldots,n-1\}\), and each integer interval \([k/n,(k+1)/n)\) has probability \(1/n\).

This is also the logic behind rng.integers(low, high).

def sample_randint_numpy(low, high, size, rng: np.random.Generator | None = None):
    """Sample from randint(low, high) using only NumPy primitives."""
    low, high = validate_low_high(low, high)
    n = high - low

    if rng is None:
        rng = np.random.default_rng()

    u = rng.random(size=size)  # Uniform on [0,1)
    return (low + np.floor(n * u)).astype(int)


low, high = 2, 7
n = 20_000
x = sample_randint_numpy(low, high, size=n, rng=rng)

values = np.arange(low, high)
emp_pmf = np.array([(x == v).mean() for v in values])
theo_pmf = np.full_like(emp_pmf, 1.0 / (high - low), dtype=float)

print("values:", values)
print("empirical pmf:", emp_pmf)
print("theoretical pmf:", theo_pmf)

values: [2 3 4 5 6]
empirical pmf: [0.1956 0.2067 0.1962 0.2079 0.1935]
theoretical pmf: [0.2 0.2 0.2 0.2 0.2]

8) Visualization

We’ll visualize:

• the PMF for several (low, high) choices

• the CDF (step function)

• Monte Carlo samples: the empirical PMF compared to the theoretical PMF

# PMF for several parameter choices
params_list = [(0, 6), (2, 8), (0, 11)]  # (low, high)

x_grid = np.arange(-1, 12)

fig_pmf = go.Figure()
for low, high in params_list:
    fig_pmf.add_trace(
        go.Bar(
            name=f"low={low}, high={high}",
            x=x_grid,
            y=randint_pmf(x_grid, low, high),
        )
    )

fig_pmf.update_layout(
    title="randint PMF for different (low, high)",
    xaxis_title="x",
    yaxis_title="P(X = x)",
    barmode="group",
)
fig_pmf.show()


# CDF for the same parameters
x_cont = np.linspace(-1.0, 12.0, 700)

fig_cdf = go.Figure()
for low, high in params_list:
    fig_cdf.add_trace(
        go.Scatter(
            name=f"low={low}, high={high}",
            x=x_cont,
            y=randint_cdf(x_cont, low, high),
            mode="lines",
            line_shape="hv",
        )
    )

fig_cdf.update_layout(
    title="randint CDF (step function)",
    xaxis_title="x",
    yaxis_title="F(x)",
)
fig_cdf.show()


# Monte Carlo: empirical PMF vs theoretical PMF
low, high = 2, 7
n = 50_000
x = sample_randint_numpy(low, high, size=n, rng=rng)

values = np.arange(low, high)
emp = np.array([(x == v).mean() for v in values])
theo = np.full_like(emp, 1.0 / (high - low), dtype=float)

fig_mc = go.Figure()
fig_mc.add_trace(go.Bar(name="empirical", x=values, y=emp))
fig_mc.add_trace(go.Bar(name="theoretical", x=values, y=theo))
fig_mc.update_layout(
    title=f"Monte Carlo check (n={n:,}) for randint(low={low}, high={high})",
    xaxis_title="x",
    yaxis_title="probability",
    barmode="group",
)
fig_mc.show()

9) SciPy Integration

SciPy provides this distribution as scipy.stats.randint.

• PMF / CDF: stats.randint.pmf, stats.randint.cdf

• Sampling: stats.randint.rvs or frozen rv = stats.randint(low, high) then rv.rvs(...)

• Entropy and moments: rv.entropy(), rv.stats(moments="mvsk")

Fitting

randint has integer shape parameters with unbounded domains, so scipy.stats.fit requires finite bounds. For this distribution, you often have an analytic MLE: $\( \hat{\ell}=\min_i x_i,\qquad \hat{h}=\max_i x_i+1. \)$

randint = stats.randint

low, high = 2, 7
print("pmf(low..high):", randint.pmf(np.arange(low, high), low, high))
print("cdf(low-1..high):", randint.cdf([low - 1, low, high - 1, high], low, high))

rv = randint(low, high)
samples = rv.rvs(size=10, random_state=rng)
print("rvs:", samples)

# Analytic MLE from a sample
data = randint.rvs(5, 13, size=2_000, random_state=rng)
low_hat, high_hat = randint_mle(data)
print()
print("Analytic MLE:")
print("  low_hat =", low_hat)
print("  high_hat=", high_hat)

# scipy.stats.fit requires finite bounds (and returns float-valued parameters)
fit_res = stats.fit(
    randint,
    data,
    bounds={
        "low": (low_hat - 5, low_hat),
        "high": (high_hat, high_hat + 5),
        "loc": (0, 0),
    },
)

print()
print("scipy.stats.fit result:")
print(fit_res)
print("  (low, high, loc) =", (fit_res.params.low, fit_res.params.high, fit_res.params.loc))

pmf(low..high): [0.2 0.2 0.2 0.2 0.2]
cdf(low-1..high): [0.  0.2 1.  1. ]
rvs: [3 3 6 2 3 3 4 5 5 6]

Analytic MLE:
  low_hat = 5
  high_hat= 13

scipy.stats.fit result:
  params: FitParams(low=5.0, high=13.0, loc=0.0)
 success: True
 message: 'Optimization terminated successfully.'
  (low, high, loc) = (5.0, 13.0, 0.0)

10) Statistical Use Cases

Hypothesis testing (goodness-of-fit to uniformity)

If you believe outcomes should be equally likely across a fixed integer set, you can test whether observed counts match a uniform distribution using a chi-square goodness-of-fit test.

Bayesian modeling (bounded discrete prior)

A discrete uniform distribution is a natural prior over a bounded set of integer hypotheses. A classic example is the German tank problem: serial numbers are modeled as i.i.d. uniform draws from \(\{1,\ldots,N\}\) with unknown \(N\).

Generative modeling (uniform mixture weights / random indices)

In generative models, you often sample an index uniformly:

• choose a mixture component uniformly

• choose a data augmentation option uniformly

• choose a random class label for synthetic data

# --- Hypothesis testing: chi-square goodness-of-fit to a known uniform support ---
low, high = 0, 10
n = 2_000
x = sample_randint_numpy(low, high, size=n, rng=rng)

values = np.arange(low, high)
counts = np.array([(x == v).sum() for v in values])
expected = np.full_like(counts, n / (high - low), dtype=float)

chi2, p_value = stats.chisquare(f_obs=counts, f_exp=expected)
print("Chi-square test for uniformity")
print("  chi2 statistic:", chi2)
print("  p-value       :", p_value)


# --- Bayesian modeling: German tank problem (posterior over N) ---
# Model: serials ~ Uniform{1,2,...,N} i.i.d. (inclusive upper bound)
# Map to SciPy's randint with low=1, high=N+1.
N_true = 200
n_obs = 15
serials = stats.randint.rvs(1, N_true + 1, size=n_obs, random_state=rng)
max_serial = int(serials.max())

N_max = 600  # prior upper limit
N_grid = np.arange(max_serial, N_max + 1)

# Uniform prior over N_grid: p(N) = constant.
# Likelihood: p(data | N) = 1/N^n_obs for N >= max_serial, else 0.
log_post_unnorm = -n_obs * np.log(N_grid)
log_post = log_post_unnorm - logsumexp(log_post_unnorm)
post = np.exp(log_post)

post_mean = float(np.sum(N_grid * post))
post_map = int(N_grid[np.argmax(post)])

cdf = np.cumsum(post)
ci_low = int(N_grid[np.searchsorted(cdf, 0.025)])
ci_high = int(N_grid[np.searchsorted(cdf, 0.975)])

print()
print("German tank example")
print("  true N        :", N_true)
print("  max observed  :", max_serial)
print("  posterior mean:", round(post_mean, 2))
print("  posterior MAP :", post_map)
print("  95% credible interval:", (ci_low, ci_high))

fig = go.Figure()
fig.add_trace(go.Scatter(x=N_grid, y=post, mode="lines", name="posterior"))
fig.add_vline(x=N_true, line_dash="dash", line_color="gray", annotation_text="N_true")
fig.add_vline(x=max_serial, line_dash="dash", line_color="red", annotation_text="max serial")
fig.update_layout(
    title="Posterior over N in the German tank problem (uniform prior)",
    xaxis_title="N",
    yaxis_title="posterior probability",
)
fig.show()


# --- Generative modeling: uniform mixture component index ---
K = 3
means = np.array([-2.0, 0.0, 3.0])
sigma = 0.6
n = 5_000

z = sample_randint_numpy(0, K, size=n, rng=rng)  # component index in {0,1,...,K-1}
y = rng.normal(loc=means[z], scale=sigma)

fig = go.Figure()
fig.add_trace(go.Histogram(x=y, nbinsx=60, name="samples"))
fig.update_layout(
    title="Samples from a 3-component Gaussian mixture with uniform weights",
    xaxis_title="y",
    yaxis_title="count",
)
fig.show()

Chi-square test for uniformity
  chi2 statistic: 10.2
  p-value       : 0.3345381516340064

German tank example
  true N        : 200
  max observed  : 196
  posterior mean: 210.54
  posterior MAP : 196
  95% credible interval: (196, 254)

11) Pitfalls

• Inclusive vs exclusive upper bound:

  • SciPy/NumPy use a half-open interval: low is included and high is excluded.

  • Python’s random.randint(a, b) is inclusive on both ends.

• Invalid parameters: you must have high > low. If high == low + 1, the distribution is degenerate (always returns low).

• Fitting requires bounds: scipy.stats.fit(stats.randint, data) fails without finite bounds because low and high are unbounded a priori.

• Non-integer parameters/observations: SciPy returns nan if low/high are not integer-valued.

• MGF overflow: \(M_X(t)\) can overflow for large positive \(t\); use it for modest \(t\) and prefer characteristic functions for numerical stability.

# Demonstrate the 'fit requires bounds' pitfall
x = stats.randint.rvs(2, 7, size=500, random_state=rng)

try:
    stats.fit(stats.randint, x)
except Exception as e:
    print("stats.fit(stats.randint, x) failed as expected:")
    print(" ", type(e).__name__ + ":", e)

# Provide finite bounds to make it work
low_hat, high_hat = randint_mle(x)
fit_res = stats.fit(
    stats.randint,
    x,
    bounds={"low": (low_hat - 3, low_hat), "high": (high_hat, high_hat + 3), "loc": (0, 0)},
)
print()
print("fit with bounds:")
print(fit_res)

stats.fit(stats.randint, x) failed as expected:
  ValueError: The intersection of user-provided bounds for `low` and the domain of the distribution is not finite. Please provide finite bounds for shape `low` in `bounds`.

fit with bounds:
  params: FitParams(low=2.0, high=7.0, loc=0.0)
 success: True
 message: 'Optimization terminated successfully.'

12) Summary

• randint(low, high) is a discrete uniform distribution on \(\{\texttt{low},\ldots,\texttt{high}-1\}\).

• PMF: \(\mathbb{P}(X=k)=1/(h-\ell)\) for integers \(k\) in the support.

• Mean: \((\ell+h-1)/2\); Variance: \(((h-\ell)^2-1)/12\); Entropy: \(\log(h-\ell)\).

• A simple NumPy sampler is \(X=\ell+\lfloor (h-\ell)U\rfloor\) with \(U\sim\mathrm{Uniform}(0,1)\).

• The likelihood for i.i.d. data is proportional to \((h-\ell)^{-m}\) when the interval contains all observations, so the MLE is \((\min x_i,\ \max x_i+1)\).

• In SciPy, scipy.stats.fit needs finite bounds for randint because parameters are unbounded.